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ABSTRACT 

We report here on variable propagation effects in over twenty years of multi-frequency timing analysis of pulsar 
PSR B 1937+21 that determine small-scale properties of the intervening plasma as it drifts through the sight 
line. The phase structure function derived from the dispersion measure variations is in remarkable agreement 
with that expected from the Kolmogorov spectrum, with a power law index of 3.66 ±0.04, valid over an 
inferred scale range of 0.2 — 50 A.U. The observed flux variation time scale and the modulation index, along 
with their frequency dependence, are discrepant with the values expected from a Kolmogorov spectrum with 
infinitismally small inner scale cutoff, suggesting a caustic-dominated regime of interstellar optics. This implies 
an inner scale cutoff to the spectrum of ~ 1 .3 x 10 9 meters. Our timing solutions indicate a transverse velocity 
of 9 km sec" 1 with respect to the solar system barycenter, and 80 km sec" 1 with respect to the pulsar's LSR. 
We interpret the frequency dependent variations of DM as a result of the apparent angular broadening of the 
source, which is a sensitive function of frequency (oc v~ 2 ,2 ). The error introduced by this in timing this pulsar 
is ~2.2 /zs at 1 GHz. The timing error introduced by "image wandering" from the slow, nominally refractive 
scintillation effects is about 125 nanosec at 1 GHz. The error accumulated due to positional error (due to image 
wandering) in solar system barycentric corrections is about 85 nanosec at 1 GHz. 

Subject headings: ISM: general — pulsars: general — radio continuum: general — scattering — turbulence 



1. INTRODUCTION 

The dispersion measure (DM) of a pulsar probes the col- 
umn density of free electrons along the line of sight (LOS). 
Observed DM variations over time scales of several weeks 
to years sample structures in the electron plasma over length 
scales of 10 m to 10 12 m. Diffraction of pulsar signals is the 
result of scattering by structures on scales below the Fresnel 
radius, 10 8 m or so. The DM as well as the scattering mea- 
sure (SM) variability along the LOS to the Crab pulsar was 
first reported by Rankin & Isaacman (1977), who reported 
that the DM variability poorly correlated with the SM vari- 
ability. Helfand et al. (1980) inferred an upper limit for DM 
variations of a few parts in a thousand for several pulsars. In 
an earlier study of PSR B 1937+21 Cordes et al. (1990) mea- 
sured a DM change of ADM ~ 0.003 pc cm" 3 over a period 
of a thousand days. The work of Phillips & Wolszczan (1991) 
presented the variations of DM observed along the LOS to 
a few pulsars. They connected these variations to those on 
diffractive scales, and derived an electron density fluctuation 
spectrum slope of 3.85 ±0.04 over a scale range of 10 7 — 10 13 
meters. Backer et al. (1993) report on further DM variability 
and show that the amplitude of the variations known at that 
time are consistent with a scaling by the square root of DM. 
Another important investigation by Kaspi et al. (1994) stud- 
ied DM variations of the millisecond pulsars PSR B 1937+21 
and B 1855+09 over a time interval of calendar years 1984 
to 1993. In addition to establishing a secular variation in DM 
over this time interval, they also show that the underlying den- 
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FIG. 1 . — Summary of our data sample. See text for details. 



sity power spectrum has an index of 3.874±0.011, which is 
close to what we would expect if the density fluctuations are 
described by Kolmogorov spectrum. An anomalous disper- 
sion event towards the Crab pulsar was reported by Backer et 
al. (2000), where they report a DM "jump" as large as 0.1 pc 
cm" 3 . 
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Table 1 

Template fit parameters at various frequencies. 



Table 2 
Parameters of PSR B1937+21. 



v (MHz) backend 


VI' 1 
'* 1 


h 


i r t 




h 


YV i 




327 


GBPP 


8.7 











186.9 


12.0 


0.51 


430 


ABPP 


10.3 


7.0 


2.5 


0.19 


186.8 


12.4 


0.53 


610 


GBPP 


9.4 


7.1 


2.9 


0.34 


187.3 


10.8 


0.60 


863 


EBPP 


8.9 


7.7 


5.2 


0.38 


187.7 


10.9 


0.55 


1000 


GBPP 


9.2 


8.6 


3.9 


0.56 


187.9 


11.3 


0.54 


1419 


ABPP 


8.5 


8.5 


3.7 


0.37 


187.5 


10.1 


0.45 


1689 


EBPP 


9.1 


9.2 


3.9 


0.37 


187.4 


10.5 


0.40 


2200 


GBPP 


9.4 


9.8 


3.2 


0.29 


187.6 


10.4 


0.36 


2379 


ABPP 


10.0 


9.8 


3.0 


0.29 


187.1 


10.8 


0.33 



Parameter 



value 



In this work, we present results of various long term moni- 
toring programs on PSR B 1937+21. Our data, which includes 
that of Kaspi, Taylor & Ryba (1994), spans calendar years 
from 1983 to 2004. These data sets have been taken with 
five different telescopes, the NRAO 1 Green Bank 42-m (140- 
ft) and 26-m (85-ft) telescopes, the NAIC 2 Arecibo telescope 
and the Nancay telescope, at frequency bands of 327, 610, 
800, 1400 and 2200 MHz. After giving the details of our ob- 
servations in 5j2] we describe our analysis methods in |3] This 
is followed in 5@]by a discussion on the distribution of scatter- 
ing material along the LOS. As we describe, the knowledge 
of temporal and angular broadening of the source, proper mo- 
tion, and scintillation based velocity estimates enables us to at 
least qualitatively study the distribution of scattering matter as 
well as properties of its wave-number spectrum. 

We have measured some of the basic refractive scintillation 
parameters from our observations, and these are discussed 
in 13 The frequency dependence of the refractive scintilla- 
tion time scale and the modulation index indicate a caustic- 
dominated regime that results from a large inner scale in the 
spectrum. 

We have detected DM variations as a function of time and 
frequency. We determine the phase structure function of the 
medium with the knowledge of the time dependent DM vari- 
ations, which is consistent with a Kolmogorov distribution of 
density fluctuations between scale sizes of about 1 to 100 A.U. 
These are summarized in ^6]and iJT] 

PSR B 1937+21 is known for its short term timing stabil- 
ity. However, the achievable long term timing accuracy is 
suspected to be seriously limited by the interstellar scatter- 
ing properties. With our sensitive measurements, we are in a 
position to quantify these errors. In |8] we describe in detail 
various sources of these errors and quantify them. 

2. OBSERVATIONS 

We have used five different primary data sets for this anal- 
ysis. The first set is the 1984-1992 Arecibo pulse timing 
and dispersion measurements obtained by Kaspi et al. (1994; 
hereafter KTR94). Their observations were performed with 
their Mark II backend (Rawley 1986; Rawley, Taylor & Davis 
1988) and later their Mark III backend (Stinebring et al. 1992) 
at two different radio frequency bands, 1420 MHz and 2200 
MHz. Their analysis methods are described in KTR94. 

' The National Radio Astronomy Observatory (NRAO) is owned and op- 
erated by Associated Universities, Inc under contract with the US National 
Science Foundation. 

2 The National Astronomy and Ionosphere Center is operated by Cornell 
University under contract with the US National Science Foundation. 
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The second data set is from 800-MHz and 1400-MHz ob- 
servations at the NRAO 140-ft telescope in Green Bank, WV. 
The Spectral Processor backend, a hardware FFT device, was 
used. Details about the observations and analysis are con- 
tained in an earlier report on dispersion measure variability 
(Backer et al. 1993). 

The third data set consists of observations at 327 MHz and 
610 MHz using the 26m (85-ft) pulsar monitoring telescope 
at NRAO's Green Bank, WV site. Room temperature (un- 
cooled) receivers at the two bands are mounted off-axis. At 
327 MHz the total bandwidth used was 5.5 MHz, and 16 MHz 
was used at 610 MHz. The two orthogonally polarized sig- 
nals were split into 32 frequency channels in a hybrid ana- 
log/digital filter bank in the GBPP (Green Bank-Berkeley 
Pulsar Processor). Dispersion effects were removed in the 
GBPP in real-time with a coherent (voltage) deconvolution 
algorithm. At the end of the real-time processing folded pulse 
profiles were recorded for each frequency channel and polar- 
ization. Further details of the backend and analysis can be 
found in Backer, Wong & Valanju (2000). PSR B1937+21 
was observed for about two hours per day starting in mid- 
1995. 

The fourth data set comes from a bi-monthly precision tim- 
ing program that includes B 1937+21 at the Arecibo Obser- 
vatory which we started in 1999 after the telescope upgrade. 
Signals at 1420 MHz and 2200 MHz were recorded using the 
ABPP backend (Arecibo-Berkeley Pulsar Processor), which 
is identical to the GBPP. Our typical observing sessions at 
1420 MHz and 2200 MHz had bandwidths of 45 MHz and 56 
MHz, respectively, and integration times of approximately 10 
minutes per session. 

The fifth data set is from a pulsar timing program that has 
been ongoing since 1989 October with the large decimetric 
radio telescope located at Nancay, France. The Nancay tele- 
scope has a surface area of 7000 m 2 , which provides a tele- 
scope gain of 1.6 K Jy _I . Observations are performed with 
dual-linear feeds at frequencies 1280, 1680 and 1700 MHz. 
Then the signal is dedispersed by using a swept frequency os- 
cillator (at 80 MHz) in the receiver IF chain. The pulse spectra 
are produced by a digital autocorrelator with a frequency res- 
olution of 6.25 kHz. Cognard et al. (1995) describe in detail 
the backend setup and the analysis procedure. 

A small amount of additional data from the Effelsberg tele- 
scope was used in our profile analysis. At Effelsberg the 



Timing PSR B 1937+21 - Interstellar Plasma Weather 



3 



EBPP backend, a copy of the GBPP/ABPP, was used. 

3. BASIC ANALYSIS 

We first present several results from the analysis of these 
data sets: a description of the frequency-dependent profile 
template used for timing; spin and astrometric timing param- 
eters from high frequency data; pulse broadening, flux den- 
sities and dispersion measure as functions of time. In ^4]we 
proceed to interpret these results and return to finer details 
regarding dispersion measure variations in 

Our basic data set consists of average pulse profiles ob- 
tained approximately every 5 minutes in each of the radio 
frequency bands - 327, 610, 800, 1420 and 2200 MHz. Fig- 
ure [0 provides a graphical summary of observation epochs 
vs date. For data sets corresponding to all frequencies ex- 
cept 327 MHz, Times of Arrival (TOAs) were computed by 
cross correlating these average profiles with a template pro- 
file. The template profile at a given frequency was made by 
using multiple Gaussian fits to very high signal to noise ra- 
tio average profiles at that frequency; the interactive program 
bfit, which is based on M. Kramer's original program fit was 
used. These fit parameters are listed in Table 1 . Col. 1 in the 
Table gives the radio frequency and the backend name is in 
col. 2. Col. 3 gives the width of component 1 (w\ \ its loca- 
tion is taken to be degrees and its amplitude is set to 1.0); 
cols. 4-6 and cols. 7-9 give the location (/), width (w) and am- 
plitude (h) values for components 2 and 3, respectively. The 
location and width are given in units of longitudinal degrees, 
where 360° indicates one full rotation cycle. The results of 
this analysis can be compared with that of Foster et al. (1991) 
which are given on the line at 1000 MHz 3 . There is reason- 
able agreement for all values except h-2 which must have been 
erroneously entered in Table 4 of Foster et al. In our analysis 
templates corresponding to arbitrary frequencies are produced 
by spline-interpolation of the component parameters. 

We used the Arecibo (1420 MHz and 2200 MHz) TOAs, 
and the GBT 140-ft (800 MHz and 1420 MHz) TOAs to fit 
for pulsar spin (rotation frequency (/), first time derivative 
(/), and second time derivative (/)) and astrometric (position 
(RAJ, DECJ), proper motion (PMRA, fj, a , along right acen- 
sion, and PMDEC, fi$, along declination) parameters. All 
TOAs were referred to the UTC time scale kept by the Na- 
tional Institute of Standards and Technology (NIST) via GPS 
satellite comparison. We removed the effects of variable dis- 
persion from this fitting procedure with weekly estimation of 
DMs and subsequent extrapolation of the dual frequency data 
to infinite frequency prior to parameter estimation. The na- 
ture of achromatic timing noise makes it particularly diffi- 
cult to determine a precise timing model. As one adds ad- 
ditional higher derivatives of rotation frequency (e.g., a third 
derivative), the best fit parameters change by amounts much 
larger than the nominal errors reported by the package that 
we used, TEMPO. The results are listed in Table 2. The er- 
rors presented in the table incorporates the range of variation 
of each parameter, as additional derivative terms are included. 
In comparison to Kaspi, Taylor & Ryba (1994), the derived 
proper motion values are marginally different. We attribute 
this difference to the variable influence of timing noise. An 
important point that needs to be stressed here is that there is 
no reason for us to assume that the higher derivative terms 
of rotation period (e.g., / or higher) has anything to do with 
the radiative braking index. They are most likely dominated 

3 The widths w\ and W3 are inverted in Table 4 of Foster et al. 
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FIG. 2. — Measured temporal pulse broadening timescale (r sc ) as a function 
of time at 327 MHz. 



by some intrinsic instabilities of the star itself, or some other 
perturbation on the star. 

Extension of dispersion measurement to 327 MHz requires 
removal of the time-variable broadening of the intrinsic pulse 
profile owing to multipath propagation in the interstellar 
medium. We deconvolved the effect of interstellar scattering 
following precepts first introduced by Rankin et al. (1970). 
We assume that the interstellar temporal broadening is quan- 
tified in terms of convolution of a Gaussian function and a 
truncated exponential function. If there is only one scatter- 
ing screen along the LOS, the assumption of a truncated ex- 
ponential function will suffice to represent the scatter broad- 
ening. However, since the scattering may arise from mate- 
rial distributed all along the LOS, a more realistic represen- 
tation is approximated by a truncated exponential function 
"smoothed" (convolved) with a Gaussian function. The in- 
trinsic pulse profile was estimated by extrapolation of param- 
eters from the higher frequency profiles. In the deconvolution 
procedure, we minimized the normalized \ 2 value by vary- 
ing the width of the Gaussian w g and the decay time scale of 
the truncated exponential r e , while keeping the intrinsic pulse 
profile fixed. The pulse scatter broadening is quantified as 
i~sc = (w 2 g + T^) l l 2 . We repeated this for average profiles ob- 
tained at every epoch to obtain the r sc measurement. In our 
fits, the average value of w g came to about 74 /isec, whereas 
the corresponding value for r e was about 85 //sec. The mea- 
surement of r sc versus time at 327 MHz is plotted in the Figure 
13 This quantity has a mean value of 120 /.is, an RMS varia- 
tion of 20 /is, and a fluctuation timescale of ~ 60 days. We 
explain these variations as the result of refractive modulation 
of this inherently diffractive parameter in discussion below. 
The estimated RMS variation at the next higher frequency in 
our data set, 610 MHz, is ^2.5 /is, using a frequency depen- 
dence of t sc ex v~ 4A . This is too small to allow fitting at this 
frequency band. 

In the strong scintillation regime, time dependent variations 
in the observed flux occur in two distinct regimes — diffrac- 
tive and refractive. The diffractive effects are dominated by 
structures smaller than the Fresnel scale, and appear on short 
time scales and over narrow bandwidths. In our observations 
diffractive modulations are strongly suppressed. On the other 



4 



R. Ramachandran et al. 



hand refractive effects occur on days time scales and are cor- 
related over wide bandwidths. We have analyzed our best data 
sets - the densely sampled data at 327 MHz and 610 MHz 
from Green Bank and at 1410 MHz from Nancay - for flux 
density variations as a function of time. The data are pre- 
sented in Figure[3] 

In analyzing the low frequency flux data from Green Bank, 
we have not adopted a rigorous flux calibration procedure. 
While there is a pulsed calibration noise source installed in 
this system, equipment changes and the nature of the au- 
tomated observing have led to large gaps in the calibration 
record. Rather than dealing with a mix of calibrated and un- 
calibrated data, or lose a large fraction of the data, we decided 
not to apply any calibration. Instead, we normalize our data 
by assuming the system temperature is constant. In order to 
see what effect this has on our results, we did two tests. 

First, we analyzed observations of PSR B1641-45, taken 
with the same system, over a similar time range. This pul- 
sar is known to have a very long refractive timescale, T le f > 
1800 days (Kaspi & Stinebring, 1992), so it can be used as a 
flux calibrator. In our data, we find it to have a modulation 
index, m = 0.10. This immediately puts a upper limit of 10% 
on any systematic gain and/or system temperature variations. 
Since modulation adds in quadrature, and we observe modula- 
tion indeces of m ~ 0.4 for PSR B1937+21, gain fluctuations 
represent at most a small fraction of the observed modulation. 

We also considered the possibility that gain variations could 
influence our measurement of T le f. This might happen if 
they occur with a characteristic timescale longer than 1 day. 
In order to test this, we analyzed observations of the Crab 
pulsar, PSR B0531+21, again taken with the same system 
over the same time range. The refractive parameters of this 
pulsar were studied in detail by Rickett & Lyne (1990). It 
makes a good comparison since it has modulation index of 
m = 0.4 at 610 MHz, very similar to PSR B1937+21. Apply- 
ing the structure function analysis (see ^5} to this data gives 
r ref = 11 days at 610 MHz, and T ref = 63 days at 327 MHz, 
consistent with the previously published results and a scaling 
law of r, e f oc v~ 12 . 

The procedure that we have adopted to calibrate our data 
set from Nancay telescope is described in detail in Cognard 
etal. (1995). 



4. DISTRIBUTION OF SCATTERING MATERIAL ALONG THE LINE 
OF SIGHT 

Several authors have shown how the scattering parameters 
of a pulsar can be used to assess the distribution of scatter- 
ing material along the LOS (Gwinn et al. 1993; Deshpande 
& Ramachandran 1998; Cordes & Rickett 1998). This re- 
sults from the varied dependences of the scattering parame- 
ters on the fractional distance of scattering material along the 
LOS. PSR B 1937+21 is viewed through the local spiral arm 
as well as the Sagittarius arm which are both potential sites 
of strong scattering. The parameters employed in this analy- 
sis are: the temporal pulse broadening by scattering (r sc ; or 
its conjugate parameter Av, the diffractive scintillation band- 
width), the diffractive scintillation time scale (Tdiff), the angu- 
lar broadening from scattering {Oh), the proper motion of the 
pulsar (/i Q , and a distance estimate of the pulsar (D). 

Let us first compare Oh and r sc that are the result of multiple 
scattering along the LOS, and express them as (Blandford & 
Narayan 1985) 



t sc = / x(D — x) ib(x) dx 

2cD J 

9 " = ^r I 



(i) 

(2) 



In these equations, x is the coordinate along the LOS, with 
the pulsar at x = and the observer at x = D. ip(x) is the 
mean scattering rate. If the scattering material is uniformly 
distributed along the LOS, then the relation between the two 
quantities can be expressed as H = 161n2 (ct sc /D). With the 
distance to the pulsar of 3.6 kpc to the pulsar (according to 
the distance model of Cordes & Lazio 2002), and the aver- 
age pulse broadening time scale of 120 fis (from the present 
work), we obtain an estimate of the angular broadening, T , of 
12 mas. This is in modest agreement with the measured value 
of 14.6±1 .8 mas (Gwinn et al. 1993), given the uncertainty in 
the distance to the pulsar and the simple assumption that the 
scattering material is uniformly distributed along the LOS. 

Next, we formulate two approaches to estimation of the 
velocity of the LOS with respect to the scattering medium, 
and use these approaches to assess the location and extent of 
the medium. The transverse velocity of the pulsar based on 
the measured proper motion (Table 2) an assumed distance of 
D = 3.6 kpc (Cordes & Lazio 2002) is 9 km s . This value 
is the velocity of the pulsar with respect to the solar system 
barycenter. With the assumed "Flat Rotation Curve" linear 
velocity of the Galaxy of 225 km s" 1 , and the Sun's peculiar 
velocity of 15.6 km s" 1 in the Galactic coordinate direction of 
(l,b) = (48.8°, 26.3°) (Murray 1986), the transverse velocity 
of the pulsar in its LSR (V p ) is 80 km s . 

The scintillation velocity (Vi ss ), which is an estimate of the 
velocity of the diffraction pattern at the location of the Earth, 
is estimated from the decorrelation bandwidth (Av) and the 
diffractive scintillation time scale (T^s), Gupta et al. (1994) 
conclude that 



V iss = 3.85 x 10 4 



IDzAu 



1 



(1-Z) Tdiff "GHz 



km s 



(3) 



where z^ghz is the observing frequency in GHz, D is in kpc, 
Av is in MHz, and Tdiff is in seconds. The variable z gives 
the fractional distance to the scattering screen, where z = 
gives the observer's position, and z = 1 gives the pulsar's po- 
sition. The value of decorrelation bandwidth is computed by 
the relation Av = l/(27rr sc ). When the effective scattering 
screen is midway along the LOS (z = 0.5), Vi ss = V p , and when 
the screen is at the location of the pulsar (z = 1.0), Vj ss = oo. 
While doing this, an important assumption is that the pulsar 
proper motion is dominant over contributions from differen- 
tial Galactic motion, solar peculiar velocity, and the Earth's 
annual orbital modulation. In the case of PSR B 1937+21, this 
assumption is not justified. The effective scattering screen, 
which is located somewhere along the LOS, has a Galactic 
motion whose component along the LOS direction is different 
from that of the pulsar or the Sun. In order to correct for this 
effect, let us calculate the LOS velocity across the effective 
scattering screen at a fractional distance z from the observer: 



V± = 3.85 x 10 



4 y/Dz(l-z)Av 



km s 



(4) 



Tdiff t'GHz 

Then, let us assume that the scattering along the LOS can 
be adequately expressed by having a thin screen alone, at a 
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FIG. 3. — Measured flux density as a function of time. The top panel corresponds to the radio frequency of 1410 MHz with the data obtained from Nancay tele- 
scope, the middle and the bottom panel to 610 and 327 MHz with the data obtained from the Green Bank 85-ft telescope. 



distance of D s = zD from the observer. Then, Equation[2]can 
be expressed as 

r sc = ^Dz(l-z) (5) 
2c 

6£ = 41n2(l-z)Vo (6) 

Here, ip gives the mean scattering rate corresponding to the 
effective thin screen. Then, let us express independently the 
transverse velocity of the LOS across the scattering screen as 

V' ± = {\-z)V e + zV p - VdzDn) 

= V E +zD(l-V c (zDn), (7) 
where V e is Earth's orbital velocity, V p is the pulsar transverse 
velocity in its LSR, and Vq is the transverse velocity contri- 
bution from the Galactic differential motion to the screen. Ve 
gives the contribution of the Earth's motion on the LOS ve- 
locity across the screen. 

Equations |4] and give two independent estimates of the 
line of sight velocity across the effective scattering screen and 
therefore allow us to solve for the value of z given D. With 
D = 3.6 kpc (Cordes & Lazio 2002), we find z = 0.7. The LOS 
velocity is 5 1 km s" 1 . The assumed value of T^ff is 78 seconds 
at 327 MHz (scaled from 444 seconds at 1400 MHz of Cordes 
et al. 1991), and the value of Ais is 0.0013 MHz calculated 
from t sc = 120 /is. 

To summarize, the measured value of Oh = 14.6 ±1.8 mas 
and the estimated value of 9 T are consistent with each other, 



suggesting a uniformly distributed scattering medium. On the 
other hand, comparison of velocity components, V p , V± and 
Vjo S suggest a thin-screen at z ~ 2/3. As Deshpande & Ra- 
machandran (1998) demonstrate, this solution is equivalent 
to having a uniformly distributed scattering medium! There- 
fore, we conclude that the line of sight to PSR B 1937+21 can 
be described adequately by a uniformly distributed scattering 
matter. 

The Earth's orbital velocity around the Sun will modulate 
the observed scintillation speed, and therefore the diffractive 
scintillation time scale, with a one-year time scale. The am- 
plitude of this modulation will depend on the effective z of the 
diffracting material, and so monitoring could provide an esti- 
mate of the effective screen location. If the effective screen 
is close to the Earth, then the modulation is strong, and if it 
is located close to the pulsar, then it is negligible. Figure [5] 
demonstrates this effect. The ordinate and abscissa give the 
LOS velocity across the effective scattering screen along the 
galactic longitude and latitude, respectively. For an assumed 
distance of 3.6 kpc, the straight line shows the expected cen- 
troid velocity of V±'. The left most end of the line (origin of 
the plot) corresponds to z = 0, and the right most end corre- 
sponds to z = 1. The annual modulation of V±', shown as the 
two ellipses, correspond to z=0.5 and z = 2/3. We have no 
way of identifying this annual modulation in our data, as we 
are insensitive to diffractive effects in our data set. 
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FIG. 4. — Structure function of normalized flux variations at 1410 MHz 
(Top), 610 MHz (middle), and 327 MHz (bottom). The 1410 MHz data was 
obtained from Nancaytelescope. The saturation value of the structure func- 
tion at larger lag values indicates the observed modulation index. Error bars 
are shown on only a few points, to preserve clarity. See text for details. 



Another measurement that could help us is the direct mea- 
surement of distance to this object by parallax measurements. 
Chatterjee et al. (2005, private communication), from their 
preliminary Very Long Baseline Array (VLBA) based paral- 
lax measurements, report that the distance to PSR B 1937+21 
is 2.3^ 5 kpc, if they force the proper motion value to be the 
same as that of our timing based measurements (Table 2). In 
the coming year, accuracy of their measurements will improve 
with further sensitive observations. 

5. REFRACTIVE SCINTILLATION 

5.1. Parameter estimation 

We determine refractive scintillation parameters from the 
data presented in Figure [3] following the structure function 
approach in previous studies (Stinebring et al 2000; Kaspi & 
Stinebring; Rickett & Lyne 1990). We define the structure 
function Dp for flux time series F(t) as 



where St is a time delay. Since our flux measurments occur 
at discrete and unevenly spaced time interals, we compute the 
flux difference for all possible lags, then average results into 
logarithmically spaced bins. 

The flux structure function typically has a form described 
by Kaspi & Stinebring et al. (1992) - a flat, noise dominated 
section at small lags, then a power-law increase which finally 
saturates at a value D s at large lags. In practice, the satura- 
tion regime may have large ripples in it, an effect of the finite 
length of any data set. In addition, the measured flux structure 
function is offset from the "true" flux structure function due 
to a contribution from uncorrelated measurement errors. At 
327 MHz and 610 MHz, we estimate this noise term from the 
short-lag (noise regime) values. At 1410 MHz (fromNancay), 
we use the individual flux error bars to get the noise level. Af- 
ter subtracting the noise value, we fit the result to a function 
of the form 



D F (St) ■■ 



D s (5t/T s ) a , [0<8t<T s ] 
D s , [St > T s ] 



(9) 



D F (6t) = 



([F(t)-F(t + St)] 2 
(F(t)) 2 



(8) 



In this fit, the power law slope a, the saturation timescale T s , 
and the saturation value D s are all free parameters. The flux 
structure function data and fits are shown in Figure |4] 

As shown in Rickett & Lyne (1990), the refractive param- 
eters can be measured from the flux structure function using 
the following relationships: The modulation index m is given 
by m = y/D s /2, and the refractive scintillation timescale T K { 
is given by Df{T k {) = D s /2. All the measured parameters, 
including those measured by earlier investigators are summa- 
rized in Table 3. 

Based on a propagation model through a simple power- 
law density fluctuation spectrum, we expect to see refrac- 
tive variations in the flux measurements on a timescale T K f ~ 
0.59fjD/V±, where V± is the line of sight velocity across the 
effective scattering screen. For the argument sake, if we as- 
sume an effective scattering screen at z = 0.5, then V± ~ 40 
km sec" 1 . With Oh = 14.6 mas, the expected refractive scintil- 
lation time scale is ^3 years at 327 MHz. This is more than an 
order of magnitude in excess of the measured value. Further- 
more, if the density fluctuations in the medium are distributed 
according to the Kolmogorov power law distribution, then the 
expected frequency scaling law is 7; e f oc A 2,2 . Our measured 
values indicate a significantly different scaling. Although it is 
consistent with T re f oc A 2 2 between 610 and 1420 MHz, it is 
not so between 327 and 610 MHz, where it is consistent with 
being directly proportional to A. Our observed modulation in- 
dex (to) values are also considerably larger than predicted, and 
show a "flatter" wavelength dependence, as listed in Table 3. 
We will address this issue in detail in the following section. 

5.2. Nature of the spectrum - inner scale cutoff 

The t hree disagreements with a simple model summarized 
in 35. 1 I f ore e us to explore a few aspects of the electron density 
power spectrum that may possibly explain what we observe. 
The effects of caustics on the observed scintillations have 
been explored by several earlier investigators, most notably 
Goodman et al. (1987) and Blandford & Narayan (1985). In 
particular, if the power law scale distribution in the medium 
is truncated at an inner scale that is considerably larger than 
the diffractive scale, as they show, the observed flux variations 
are dominated by caustics. This is of great interest to us, as 
this seems to explain all the discrepancies that we note in our 
observed refractive parameters. For instance, as Goodman et 
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FIG. 5. — Estimated line of sight transverse velocity across the effective 
scattering screen at a distance of zD from the Sun. The velocity is resolved 
into the components along galactic longitude (V/) and latitude (Vj,). We have 
assumed a distance of 3.6 kpc to the pulsar from the Sun. Any point on the 
line indicates a combination of Vj and V/, corresponding to a value of z, with 
the left most end for z = and the right most end for z = 1 ■ The line itself does 
not include the Earth's orbital velocity contribution. The annual modulation 
due to Earth's motion in its orbit is shown as the ellipses. The two ellipses 
correspond to the scenarios of z = 0.5 and z = 2/3. 



objects have a strong thin-screen scatterer somewhere along 
the LOS. This is either a supernova remnant (or a plerion) 
like in the case of Vela and Crab pulsars, a HII region as in 
the case of B 1942-03 and B 1642-03, or a Wolf-Revet star 
as in the case of B 1933+16 (see Prentice & ter Haar 1969; 
Smith 1968). Although our measurements show that pulsar 
PSR B 1937+21 is consistent with the characteristics of the 
super-Kolmogorov group, as we describe in ^ we find no 
compelling evidence for the presence of any dominant scat- 
terer somewhere along the LOS. 

To summarize, while some investigators have reported 
agreement of the measured refractive properties with the theo- 
retical expectations from a Kolmogorov spectrum with an in- 
finitismally small inner scale, there are a considerable number 
of cases where the observed properties are significantly dif- 
ferent from that predicted by the simple Kolomogorov spec- 
trum. These other cases can be explained by invoking spec- 
trum with a large inner scale cutoff, including the case where 
the cutoff approaches the Fresnel radius and leads to a caustic- 
dominated regime. From Gupta etal. (1993) and Stinebring et 
al. (2000), the modulation index can be specified as a function 
of the inner cutoff scale as 

With the known value of at 327 MHz of 1.33 kHz, the 
distance to the pulsar of 3.6 kpc, and the observed modulation 
index of 0.39, the inner scale cutoff, r,, comes to 1.3 x 10 9 
meters. 



al. (1987) show that if the inner scale cutoff is a consider- 
able fraction of the Fresnel scale, then the observed fluctu- 
ation spectrum of flux is dominated by fluctuation frequen- 
cies that are lower than the diffractive frequencies, but signif- 
icantly higher than that expected from refractive scintillation. 
This is what we observe. Moreover, as they note, the observed 
wavelength dependence of the refractive time scale, as well as 
the modulation index is expected to be "shallower" than the 
expected values of A 22 and A -0 57 , respectively. 

A "shallow" frequency dependence of the modulation in- 
dex has been reported by others (Coles et al. 1987; Kaspi & 
Stinebring 1992; Gupta et al. 1993; Stinebring et al. 2000). 
While Kaspi & Stinebring (1992) find that the observed re- 
fractive quantities matched well with the predicted values for 
five objects, three other objects, especially PSR B0833^1-5, 
has a significantly shorter measured T K f and greater modula- 
tion index than expected. This is very similar to our situation 
here with PSR B1937+21. 

Stinebring et al. (2000) concluded that the 21 objects that 
they analysed fell into two groups. The first group followed 
the frequency dependence predicted by a Kolmogorov spec- 
trum with the inner cutoff scale far less than the diffractive 
scales (Kolmogorov -consistent group). The second group, 
which is the super-Kolmogorov group, is consistent with a 
Kolmogorov spectrum with an inner scale cutoff at ~ 10 8 
meters. The observed modulation indices were consistently 
greater than that of the Kolmogorov predictions, as we have 
seen in our measurements of PSR B 1937+21. This group in- 
cludes pulsars like PSRs B0833-45 (Vela), B0531+21 (Crab), 
B0835^1, B191 1-04 and B1933+16. An important physical 
property that binds them all is that, excepting one object, all 



6. DM VARIATIONS 

We turn now to the dispersion measure variations presented 
in Figure|6]that sample density variations on transverse scales 
much larger than those involved with diffractive and refractive 
effects. The most striking feature in Figure[6]is the large sec- 
ular decline from 71.040 pc cm" 3 in 1985 to 71.033 pc cm -3 
in 1991 and then to 71.022 pc cm" 3 by late 2004. These long- 
term secular variations are many times greater than the RMS 
fluctuations of ~ 10~ 3 pc cm" 3 on short time scales. An im- 
portant question that arises is whether these variations are the 
result of a spectrum of electron-density turbulence, or whether 
there might be a contribution from the smooth gradient of a 
cloud, or clouds along the LOS. We look at this question from 
two angles. First we present a phase structure function anal- 
ysis of the dispersion measure data and estimate a power-law 
index of the electron density spectrum. Then we estimate the 
probability that such a spectrum would produce a 22-y real- 
ization that was so strongly dominated by the large, mono- 
tonic changes mentioned above. 

We write the power spectrum of electron-density fluctua- 
tions as 

P(q) = C 2 n q-P, [q <q<qi] (11) 

where j3 is the power law index, q and q, are the spatial fre- 
quencies corresponding to the outer and the inner boundary 
scale, within which this power law description is valid. C 2 is 
the amplitude, or strength, of the fluctuations. A quantity that 
is closely related to the density spectrum which can be quan- 
tified by observable variables is the phase structure function, 
D^ib), with b = 2ir/q. This is defined as the mean square ge- 
ometric phase between two straight line paths to the observer, 
with a separating distance of b between them in the plane nor- 
mal to the observer's sight line. The phase structure function 
and the density power spectrum are related by a transform 
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Table 3 

Measured and expected parameters. 
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?diff 


r rc f 
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V 


(flS) 
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(s) 


observed 


expected 


observed 


expected 


(MHz) 


120+ 


14.6* 




73 days 


3/ 


0.33 


0.1 4 C 


327 


38 e 




10() c 










430 








43.9 days 


6 mon g 


0.39 


0.2 


610 






260" 


3 days'"' 


45 daysS 


0.45 


0.33 c 


1400 


0.1 7 e 




444 e 










1400 



t Has a time dependent RMS variation of 20/.ts 
"Cordes et al. (1986) 
*Gwinnetal. (1993) 

c Romani et al. (1986); Kaspi & Stinebring (1992) 

rf Lestrade et al. (1998) give the value as 13 days 

f Cordes et al. 1990 

^Calculated with 7/ rcf ~ 8 H D/2V ± 

g Extrapolated with T Icl oc A 22 



(Rickett 1990; Armstrong, Rickett, Spangler 1995), 



poo />oo 

D^b) = / 87r 2 A 2 r e 2 dz' / q [1 
Jo Jo 



-Jo(bqz'/z)] dq 



xP(q = 0) (12) 

Here, r e is the classical electron radius (2.82 x 10~ 15 meters), 
J is the Bessel function. Under the conditions that we have 
assumed, D<p{b) is also a power law (Rickett 1990; Arm- 
strong, Rickett & Spangler 1995), given by 



D*(b) = — 



'coh 



0-2 



(13) 



where b co ^ is the coherence spatial scale. Dispersion measure 
can be written as 



DM = 2.410 x 10" 



'1"2 



It J pccm ' 



(14) 



where Atfi is the difference in the arrival phases (<fc- </>i) of 
the pulse at the two barycentric radio frequencies (Hz) v\ and 
V2, with / being the barycentric rotation frequency (Hz) of 
the pulsar. With this linear relation between DM and geomet- 
ric phase difference, then structure function can be written as 
(KTR94) 

D^{bo) = 



16 pc cm 3 



(15) 



2.410x10" 

x([DM(b + b )-DM(b)] 2 ) 

Here, the angular brackets indicate ensemble averaging. The 
transformation between the spatial coordinate b (and the spa- 
tial delay b a ) and the time coordinate t (or time delay r) is 
simply given by b = V±t, where V± is the transverse velocity 
of the LOS across the effective scattering screen. 

With the understanding that any difference in DM that we 
compute for a time baseline from Figure [6] corresponds to a 
point in the phase structure function, we can derive the phase 
structure function on the basis of Equation^] This is given 
in Figure 

There are several important points in Figure The solid 
line gives the best fit line for the data in the time interval of 5 
days to 2000 days. The derived values of the intercept and the 
power law index (/3) are, 

= 4.46 ±0.09 



intercept = 



/3= 3.66 ±0.04 



(16) 



The value of /3 is remarkably close to the value expected from 
a Kolmogorov power law distribution (/3 = 1 1 /3). We are us- 
ing the terminology "intercept" only to indicate the value of 
log[D0(r)] when log[timelag(days)] is zero. Here, a caution- 
ary remark is warranted. Given the finite time span of our data 
set, and the fact that the low spatial frequencies dominate the 
long term variations in DM, we do not have a stationary sam- 
ple of noise spectrum. We have estimated the error in each 
bin of the structure function as 

- °" D 

° S ~ y/WC 

where od is the root mean square deviation with respect to 
the mean phase structure function value in a bin, D<f,(f), and 
Ni is the number of "independent" samples in the bin. This 
is estimated as the smaller of (T/t) and the actual number of 
samples that have gone into the estimation of D^if). Here, T 
is the time span of the data. 

By assuming that the transverse speed of the sightline 
across the effective scattering screen is ^40 km sec -1 (half of 
pulsar's velocity in its LSR), we can translate the delay range 
between which this slope is valid, to 0.2 to 50 A.U. 

The time delay value corresponding to the phase structure 
function value of unity is, by definition, the coherent diffrac- 
tive time scale (Tdiff) at the corresponding radio frequency, 
with the assumption that the scattering material is uniformly 
distributed along the LOS. From the fit parameters given in 
Eauationll6l this delay is 180 seconds. This should be com- 
pared with the measured Tag value of 444±28 seconds tabu- 
lated in Table 3. If we interpret the inner scale cutoff value 
of r, ; ~ 1.3 x 10 9 meters as the scale size below which the 
slope (/?) of the density fluctuation spectrum changes to a 
value greater than that given in Equation^! then the fact that 
the measured Tas value of 444 seconds being significantly 
greater than 180 seconds is understandable. In the limiting 
case, where the slope of the density irregularity power spec- 
trum changes to the critical value of (3 = 4 below the inner 
scale cutoff value, the expected Tub value is about 1 100 sec- 
onds. This makes it very important to measure the exact fre- 
quency dependence of the diffractive parameters like tempo- 
ral scatter broadening and diffractive scintillation time scale. 
To the best of our knowledge, Cordes et al. (1990) show the 
most complete multi-frequency measurements of the diffrac- 
tive scintillation parameters of this pulsars. Their measure- 
ments are not accurate enough to distinguish between such 
small variations in slope. 

While our analysis of DM variability suggests a Kol- 
mogorov spectrum at AU scales, we are struck by the long 
term monotonic decrease of DM and wonder if we might be 
seeing the effects of smooth gradients in large scale galac- 
tic structures that are not part of a turbulent cascade. We 
performed a Monte Carlo simulation to investigate this. In 
each realization of the simulation, we generated with a dif- 
ferent random number seed, a screen of density fluctuations. 
We assumed that the random fluctuations corresponding to a 
given spatial frequency are described by a Gaussian function, 
but the total power as a function of spatial frequencies is de- 
scribed by a single power law of index -1 1/3. Assuming that 
the screen is located at the mid point along the sight line, we 
let the pulsar drift with its transverse velocity, and measured 
the implied column density (DM) as a function of time. 

We developed a procedure similar to that of Deshpande 
(2000) to compare the observed DM(t) curve with the simu- 
lated ones. From the observed DM(t) curve, we computed the 
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FIG. 6. — Dispersion Measure as a function of time. Open triangles give the measurements of KTR94 at 1400 and 2200 MHz, open circles are from our Green 
Bank 140-ft telescope measurements at 800 and 1400 MHz, and the open diamond symbols indicate our measurements from the Arecibo Observatory, at 1420 
and 2200 MHz. All error bars indicate RMS errors. 
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FIG. 7. — The phase structure function derived with the help of Eguation ll5l from the data displayed in Figure^] Solid line represents the best fit for the data in 
the time range of 30 days to 2000 days. For translating the time delay range into a space delay, we have assumed a sightline transverse velocity of 40 km sec -1 , 
which is half of the pulsar's peculiar velocity in its LSR. We have assumed an effective screen at a fractional distance of 0.5. See text for details. 
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parameter ADM = [DM(t)—DM(t—T )], where r is the time 
delay. Our aim is to compare the distribution of this parameter 
in very short delays and very large delays. As we can see in 
Figure0 the structure function describes a well defined slope 
between the delay range of ~30 days to ~2000 days. We 
defined two delay bins, 30-60 and 1300-2000 days, within 
which we monitored the distribution function of the quantity 
ADM. From this, we could infer that the distribution at the 
bin of 1300-2000 days had a span of ~ 20er s , where a s is the 
RMS deviation of the distribution at the delay range of 30-60 
days. That is, the ADM values that we see at largest delays is 
as high as 20 times that of the typical deviations at short de- 
lays. We performed the same procedure on the simulated set 
of data to quantify the likelihood of such deviations. We simu- 
lated 1024 number phase screens. Out of these 1024 screens, 
we found that such large deviations were possible ~7% of 
the times. This is perhaps not surprising, as with such a steep 
spectrum, it is obvious that most of the power is in large scales 
(smaller spatial frequencies), and hence they tend to dominate 
our ADM measurements. We conclude that while monotonic 
changes of this magnitude are rare, it is consistent with a tur- 
bulent cascade spectrum of density fluctutations. 

7. FREQUENCY DEPENDENCE OF DM 

Dispersion depends on the column density of electrons. In 
a uniform medium radio wave propagation senses the aver- 
age density in a tube whose beam waist is set by the Fres- 
nel radius \/z(\ -z)XD. In a turbulent medium frequency- 
dependent multipath propagation can expand this tube consid- 
erably. With refraction, the center of the tube wanders from 
the geometric LOS. Indeed there may be a number of wave 
propagation tubes, each with their independent relative gain. 
The consequence is that DM and related effects will show fre- 
quency dependence: 

1 . the effective DM depends on frequency. 

2. the DM variations at lower frequencies will be much 
"smoother" than that at higher frequencies, as the ap- 
parent angular size of the source acts as a smoothing 
function on the measured DM variations. 

3. since the apparent size of the source is larger at low fre- 
quencies, some features of the ISM that are visible at 
lower frequencies may be invisible at higher frequen- 
cies! 

We can explore these effects by assuming that the timing 
residuals at 327 MHz and 610 MHz, which are relative to 
the timing model derived at higher frequencies that included 
removal of DM variations, are due to DM variations. The 
smoothing effect of scattering could be revealed by a spectral 
analysis. The slow variations were removed to pre-whiten the 
spectrum that would otherwise be severely contaminated. The 
327 MHz data was fit to a fourth order polynomial and the re- 
sult was subtracted from both data sets. The two right side 
panels give the residual DM values after subtracting the best 
fit curve from the actual DM curve. The resulting spectral 
comparison fails to have sufficient signal to clearly demon- 
strate increased smoothing at 327 MHz relative to 610 MHz. 
Higher signal-to-noise ratio is required. The DM variations 
relative to long term trends in the right-hand panels of Figure 
[8] are different. 

An important source of systematic error that can affect our 
analysis here is the effect of scattering on the derived DM as 



a function of time at a given frequency. At 327 MHz, as we 
described in |5] we perform an elaborate procedure to fit for 
the scatter broadening of the pulse profile, in order to compute 
the "true" TOA of the profile. However, we do not follow this 
procedure at 610 MHz (or any other higher frequency). The 
error due to this can be quantified easily from Figure [2] The 
temporal scatter broadening value varies by an RMS value of 
some 19.6 /is. With the wavelength dependence of r sc oc A 4,4 , 
the expected RMS variation at 610 MHz is 1.3 /.ts. The equiv- 
alent DM perturbation at 610 MHz with respect to infinite 
frequency is ~ 10~ 4 pc cm" 3 . 

8. ACHIEVABLE TIMING ACCURACY 

In this section, we will estimate quantitatively errors in- 
troduced by various scintillation related effects. For PSR 
B 1937+21, although a typical observation with highly sensi- 
tive telescopes like Arecibo telescope helps us achieve a TOA 
accuracy of a few tens of nanoseconds, the ultimate long term 
accuracy seems to be far greater than this. In general, it is a 
combination of frequency independent "intrinsic timing nose" 
from the pulsar itself, and the frequency dependent effects, 
such as what we are addressing here. With some 18 years of 
data at 800, 1400 and 2200 MHz, Lommen (2002) quantifies 
the timing timing residual, after fitting for position, proper 
motion, rotation frequency and its time derivative (see also 
Kaspi et al. 1994). A large fraction of the left over residuals 
is presumably the intrinsic timing noise. As we have men- 
tioned before, we have absorbed a good part of this by fitting 
for the second time derivative of the rotation frequency, / (see 
Table 2). In this section, our aim is to quantify possible timing 
errors from various "chromatic" effects related to interstellar 
scintillation. 

8. 1 . Fluctuation of apparent angular size 

The temporal variability of pulse broadening, r sc , (as shown 
in Figure[2) means that even the apparent angular broadening 
of the source, Oh, is also changing as a function of time. Since 
Tsc oc H , with the RMS variation in t sc of 19.6 /is at the ra- 
dio frequency of 327 MHz, the corresponding variation in Oh 
comes out to be ~8% of the mean value. This change occurs 
with typical time scales of ~67 days, which is the time scale 
with which r sc changes. Since we have only one epoch of Oh 
measurement, we have no way of observationally verifying 
the mean value or the time scale of its variation. 

8.2. Image wandering and the associated timing error 

Due to non-diffractive scintillation that "steers" the direc- 
tion of rays ("refractive focussing"), the position of the pulsar 
is expected to change as a function of time. This is an impor- 
tant and significant effect, as it introduces a TOA residual as 
a function of time, depending on the instantaneous position 
of source on the sky. Several authors have investigated this 
effect in the past (Cordes et al. 1986; Romani et al. 1986; 
Rickett & Coles 1988; Fey & Mutel 93; Lazio & Fey 2001). 
For a Kolmogorov spectrum of irregularities (J3 = 11/3) with 
infinitismally small inner scale cutoff, Cordes et al. (1986) 
predict the value of RMS image wandering as 

(60 2 )^ 2 = 0.18 mas (^) "\f (17) 
V A cm / 

For an assumed distance to PSR B1937+21 of 3.6 kpc, this 
comes to 2 mas at 327 MHz (wavelength, A = 92 cm). The 
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FIG. 8. — Dispersion Measure variations at 327 and 610 MHz. The two left hand side panels give DM as a function of time. The 327 MHz best fit line, 
produced by a fourth order polynomial fit, is given as a solid line in both 327 and 610 MHz plots. The residual DM, which is the difference between the actual 
DM and the best fit line, is given in the right hand side panels. See text for details. 



value of 2 mas is still significantly less than the apparent an- 
gular size of the source, 14.6 mas, measured by Gwinn et al. 
(1993). However, for a spectrum with a steeper slope or with 
a significantly larger inner scale cutoff (as in our case), the 
value of {SO 2 ) 1 ! 2 is expected to be much larger, perhaps com- 
parable to the value of 6h- 

In order to estimate the timing errors introduced by this 
image wandering, we need an estimate of scattering measure 
(SM) and C 2 along the LOS to this pulsar. Following Cordes 
etal. (1991), 



SM: 



-[ cl 

Jo 



(x) dx 



128 mas 



5/3 



11/3 
7 GHz 



:292 



D 



kpc 



5/6 



11/3 
GHz ' 



(18) 



Here, r sc is specified in seconds. SM is specified in units of 
kpc m -20 / 3 . Assuming a distance of 3.6 kpc, r sc = 120 /is, 
^=0.327 GHz, the value of SM comes to - 8.8 x 10" 4 kpc 
m -20/3 Assuming that the scattering material is uniformly 
distributed along the LOS, C 2 ~ 2.4 x 10" 4 m -20 / 3 . 

Then, for a Kolmogorov spectrum, the RMS timing residual 
due to the image wandering can be written as (Cordes et al. 



1986) 



<7 5te =26.5 ns ^ 49 / 15 D 2 / 3 ( 



C 2 



I 10- 4 m- 20 / 3 



4/5 



(19) 



With the computed value of C 2 and a distance of 3.6 kpc, this 
amounts to 4.8 fis at 327 MHz. Given the frequency depen- 
dence, this effect can be minimized by timing the pulsar at 
higher frequencies. For instance, at frequencies of 1 GHz and 
2.2 GHz, this error translates to 125 and 2 nanosec, respec- 
tively. However, given the significantly large value of the in- 
ner scale cutoff, the RMS timing error that we have computed 
may well be a lower limit, and it is likely to be higher. 

Given the fact that the exact source position due to this ef- 
fect is unknown at any given time, it is very difficult to com- 
pensate for this effect. 

8.3. Positional errors in solar system barycentric corrections 

As we saw above, due to image wandering, the apparent 
position of the source wanders in the sky. This introduces yet 
another timing error. While translating the TOA at the obser- 
vatory to the solar system barycenter, we assume a position 
which is away from the actual apparent position at the time of 
observation. This introduces an error, which can be quantified 
as (Foster & Cordes 1990) 



Afbaxy = ^(rZ-n)(l-z)A6 r (\). 



(20) 
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Here, c is the velocity of light, the dot product term gives the 
projected extra path length travelled by the ray due to Earth's 
annual cycle around the Sun, and A9 r (X) is the positional er- 
ror due to image wandering. Of course, this term is a function 
of frequency, and hence the error accumulated is different at 
different frequencies. 

At 327 MHz, with an RMS image wandering angle of 2 
mas, for an object at the ecliptic plane, Afbary ~ 2 /US. For PSR 
B1937+21, this error amounts to ^0.8 /as. At frequencies 1 
GHz and 2.2 GHz, this error translates to 85 and 17 nanosec, 
respectively. 

8.4. Timing error due to DM variation as a function of 
frequency 

An important issue that arises due to the frequency depen- 
dent DM variation is the timing accuracy. Pulsars like PSR 
B 1937+21 are known for the accuracy to which one can com- 
pute the pulse TOA. Given this, one wishes to eliminate any 
error that is incurred due to systematic effects like what we 
have here. Between 327 and 610 MHz (the two curves in 
Figure[3), the typical relative fluctuation of DM that we see is 
about 5 x 1CT 4 pc cm" 3 . As we discussed before, this is purely 
due to the fact that the effective interstellar column length 
sampled at one frequency is different from that at another fre- 
quency, due to the scatter broadening of the source. At 610 
MHz, this relative DM fluctuation corresponds to some 6 /is 
at 610 MHz. That is, at 610 MHz, typically an unaccounted 
residual of 6 /is is incurred due to just effective DM errors. 
Even if the behavior of the pulse emission is extremely stable, 
at low frequencies, interstellar scattering limits our timing ca- 
pabilities. 

Due to the fact that dispersion delay goes as v~ 2 , although 
the above mentioned effect seems significant, one should be 
able to reduce it by going to higher frequencies. For instance, 
at 2.2 GHz, the DM-limited TOA error for PSR B1937+21 
will be ^0.5/is. This is not necessarily encouraging news, as 
a timing residual error of 0.5/is is large when compared to 
the accuracy that we can achieve in quantifying the TOAs (a 
few tens of ns) for this pulsar, given our observations with 
sensitive telescopes like Arecibo. 

To summarize, although one takes into account time depen- 
dent DM changes while analysing the data, in order to achieve 
high accuracy timing, it is important to correct for a frequency 
dependent DM term. This introduces another dimension of 
correction in the timing analysis. 

9. CONCLUDING REMARKS 

We have presented in this paper a summary of over twenty 
years of timing of PSR B 1937+21. These observations have 
been done over frequencies ranging from 327 MHz to 2.2 
GHz with three different telescopes. 

Given the agreement between the measured apparent angu- 
lar broadening and that estimated by the temporal broadening, 
and the measured proper motion velocity and that estimated 



by the knowledge of scintillation parameters, we conclude 
that the scattering material is uniformly distributed along the 
sightline. 

There are three significant discrepancies between the ex- 
pected values and the measured refractive parameters. These 
are, 

1 . The measured flux variation time scale is about an or- 
der of magnitude shorter than what is expected from 
the knowledge of the observed apparent angular broad- 
ening. 

2. The flux variation time scale is observed to be directly 
proportional to the wavelength, whereas it is expected 
to vary as proportional to A 2 2 (for a Kolmogorov spec- 
trum). 

3. The flux modulation index is observed to have a wave- 
length dependence that is much "shallower" than the 
expected value. 

These three discrepancies consistently imply that the optics 
is "caustic-dominated". This would mean that the density ir- 
regularity spectrum has a large inner scale cutoff, 1.3 x 10 9 
m. Our extrapolation of the phase structure function from the 
regime sampled by DM variations to the diffractive regime 
seems to indicate that the expected Tdiff value is considerably 
shorter than the measured value. This is in favor of the above 
conclusion. Accurate measurements of frequency dependence 
of diffractive parameters is much needed. 

In general, Millisecond pulsars are known for their timing 
stability. Potentially, we may achieve adequate accuracy in 
timing some of these pulsars to understand some of the most 
important questions related to the gravitational background 
radiation, or the internal structure of these neutron stars. How- 
ever, our analysis here shows that interstellar scattering could 
be an important and significant source of timing error. As we 
have shown, although PSR B 1937+21 is known to produce 
short term TOA errors as low as 10-20 nanosec with sensitive 
observations, the long term error is far larger than this. After 
fitting for / (which absorbes most of the achromatic timing 
noise), the best accuracy that we can achieve for this pulsar 
is 0.9 /isec at 1.4 GHz, and about 0.5 /isec at 2.2 GHz (by 
one of the authors, Andrea Lommen). It appears that almost 
all of this error can be accounted for by various effects that 
we have discussed in |8] In general, for millisecond pulsars 
with substantial DM, even if achromatic timing noise is small, 
interstellar medium may be a major source of timing noise. 

We thank M. Kramer for sharing the data from the 
Effelsberg-Berkeley Pulsar Processor (EBPP), and S. Chatter- 
jee and W. Brisken for sharing their VLBA based proper mo- 
tion and parallax results of PSR B 1937+21 prior to the pub- 
lication. This work was in part supported by the NSF grant, 
AST-9820662. 
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